clear;
clc;
tstart=1; %start time
tmax=500; %max time of simulation
time=tstart:tmax; %time vector

%nparticles = 1000;
nparticles = 1;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=1:npos;
nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(1,npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(1,npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(1,npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(1,npos);%a's for spin down at each N at current time

%sigma=zeros(1,npos); %the standard deviation of the spatial distribution
sigma=tmax/2:100:tmax;
sigma_count = 1;
T=time;
pos_counter=0;
a=0;

alpha=0.025;
avg=0;
alpha1=0;
alpha2=0;
alpha3=0;
A1=0;
A2=0;
B1=0;
B2=0;

%p_n_tmax_by_2=zeros(1,npos);
p_n_tmax = zeros(1,npos);

%for all particles
for p=1:nparticles
    %form the initial a vectors for time = 0. consider all the particles to be
    %with spin up
    a=rand;
    a_n_u_prev(nmid_rel)=a;
    a_n_d_prev(nmid_rel)=1-a;
    
    %for all times
    for t=time
        display(t);
        pos_counter = 0;
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
        
        %Set the alphas for this nrange
        alpha=zeros(1,npos)+0.025;
        avg=zeros(1,npos);
        alpha1=normrnd(avg,alpha);
        alpha2=normrnd(avg,alpha);
        alpha3=normrnd(avg,alpha);
        
        %for all possible positions for this particular time, fill up the a
        %arrays for the current time
        
        for pos_counter=nstart_rel:nend_rel
%             //a=1/(sqrt(2));
            a=rand;
            
            A1= a*exp(i*alpha3(pos_counter)) + sqrt(1-a*a) * exp(-alpha2(pos_counter)+i*alpha1(pos_counter));
            A2 = a* exp(alpha2(pos_counter)+i*alpha1(pos_counter)) + sqrt(1-a*a) * exp(-i*alpha3(pos_counter));
            B1 = sqrt(1-a*a) * exp(i*alpha3(pos_counter)) - a * exp(-alpha2(pos_counter)+i*alpha1(pos_counter));
            B2 = sqrt(1-a*a) * exp(alpha2(pos_counter)+i*alpha1(pos_counter))- a * exp(-i*alpha3(pos_counter));
            
            %hadamard with noise
            a_n_u_cur(pos_counter)=  ( A1 * a_n_d_prev(pos_counter-1) )+ ( A2 * a_n_u_prev(pos_counter-1) );
            a_n_d_cur(pos_counter)=  ( B2 * a_n_u_prev(pos_counter+1)) + ( B1 * a_n_d_prev(pos_counter+1) );
            
            
            %hadamard without noise
%             a_n_u_cur(pos_counter)=  ( sqrt(1-a*a) * a_n_d_prev(pos_counter-1) )+ ( a * a_n_u_prev(pos_counter-1) );
%             a_n_d_cur(pos_counter)=  ( sqrt(1-a*a ) * a_n_u_prev(pos_counter+1)) - ( a * a_n_d_prev(pos_counter+1) );

              %eqns from paper
%             a_n_u_cur(pos_counter) =  (1/sqrt(2))*( a_n_d_prev(pos_counter-1) +  a_n_u_prev(pos_counter-1) );
%             a_n_d_cur(pos_counter)= (1/sqrt(2))*( a_n_u_prev(pos_counter+1) - a_n_d_prev(pos_counter+1) );

            %if time is tmax/2 or tmax then calculate the p_n arrays as average
            %of a arrays 
%             if(t ==tmax)
%             if(t >= tmax/2 && mod(t,100)==0)
               p_n_tmax(pos_counter) = a_n_u_cur(pos_counter)^2 + a_n_d_cur(pos_counter)^2; 
%             end;
            
        end;
        
        %calculate the sigma 
%         if(t==tmax)
%         if(t >= tmax/2 && mod(t,100)==0)
%             sigma = sqrt( (pos^2 .* p_n_tmax) - (pos.*p_n_tmax)^2 )
            sigma(sigma_count) = sqrt( sum(pos.^2 .* p_n_tmax) - (sum(pos.*p_n_tmax)).^2 );
            sigma_count = sigma_count + 1;
%         end;
        
        %set the a arrays of prev time with values of that of the current
        %time and a arrays of current time with zero values
%         if(t <tmax)
            a_n_u_prev = a_n_u_cur;
            a_n_d_prev = a_n_d_cur;
            a_n_u_cur = zeros(1,npos);
            a_n_d_cur = zeros(1,npos);
%         end;
        
    end;
    
end;

figure(1)
plot(pos-nmid_rel,p_n_tmax);

figure(2)
plot(T,sigma);
grid on;

